﻿#pragma once

#include <cassert>
#include <deque>
#include <vector>

namespace CAGD
{

/**
 * 多项式模板类.
 */
template <class T> class TPolynomial
{
  protected:
    std::vector<T> m_coefficients; // 系数向量

  public:
    // 限定显式调用，防止 operator* 中的隐式转换
    explicit TPolynomial(int o = 0) : m_coefficients(std::vector<T>(o, T()))
    {
    }

    TPolynomial(std::initializer_list<T> l)
    {
        for (auto it = l.begin(); it != l.end();)
            m_coefficients.push_back(*it++);
    }

    TPolynomial(const TPolynomial &poly) : m_coefficients(poly.m_coefficients)
    {
    }

    void operator=(const TPolynomial &poly)
    {
        m_coefficients = poly.m_coefficients;
    }

    bool operator==(const TPolynomial &poly)
    {
        return m_coefficients == poly.m_coefficients;
    }

#define POLY_OP(func, op)                                                                                              \
    TPolynomial func(const TPolynomial &poly) const                                                                    \
    {                                                                                                                  \
        int index = 0;                                                                                                 \
        int O1 = Order();                                                                                              \
        int O2 = poly.Order();                                                                                         \
        int Maxsize = std::max(O1, O2);                                                                                \
        int Minsize = std::min(O1, O2);                                                                                \
                                                                                                                       \
        TPolynomial res(Maxsize);                                                                                      \
        while (index < Minsize)                                                                                        \
        {                                                                                                              \
            res.m_coefficients[index] = m_coefficients[index] op poly.m_coefficients[index];                           \
            index++;                                                                                                   \
        }                                                                                                              \
                                                                                                                       \
        while (index < O1)                                                                                             \
        {                                                                                                              \
            res.m_coefficients[index] = m_coefficients[index];                                                         \
            index++;                                                                                                   \
        }                                                                                                              \
                                                                                                                       \
        while (index < O2)                                                                                             \
        {                                                                                                              \
            res.m_coefficients[index] = res.m_coefficients[index] op poly.m_coefficients[index];                       \
            index++;                                                                                                   \
        }                                                                                                              \
        return res;                                                                                                    \
    }
    POLY_OP(operator+, +)
    POLY_OP(operator-, -)

#undef POLY_OP

#define POLY_OP(func, op)                                                                                              \
    TPolynomial func(T k) const                                                                                        \
    {                                                                                                                  \
        int O = Order();                                                                                               \
        TPolynomial res(O);                                                                                            \
        for (int i = 0; i < O; i++)                                                                                    \
            res.m_coefficients[i] = m_coefficients[i] op k;                                                            \
        return res;                                                                                                    \
    }

    POLY_OP(operator*, *)
    POLY_OP(operator/, /)
#undef POLY_OP

    TPolynomial operator*(const TPolynomial &poly) const
    {
        const int O1 = Order();
        const int O2 = poly.Order();
        TPolynomial res(O1 + O2 - 1);
        for (int i = 0; i < O1; i++)
            for (int j = 0; j < O2; j++)
                res.m_coefficients[i + j] += m_coefficients[i] * poly.m_coefficients[j];
        return res;
    }

    TPolynomial operator/(const TPolynomial &poly) const
    {
        int O1 = Order();
        int O2 = poly.Order();

        // 阶数不够，返回 0 多项式
        if (O1 < O2)
            return TPolynomial(0);

        TPolynomial res(O1 - O2 + 1);

        // 最高阶系数
        res.m_coefficients[O1 - O2] = m_coefficients[O1 - 1] / poly.m_coefficients[O2 - 1];

        // 获得余式
        TPolynomial tmp = *this - res * poly;

        // 减去最高项后需要降阶，然后递归计算商
        tmp.Reduce();
        tmp = tmp / poly;

        // 商式赋值
        for (int i = 0; i < tmp.Order(); ++i)
            res.m_coefficients[i] = tmp.m_coefficients[i];

        return res;
    }

    // 计算多项式的值
    T Value(T X) const
    {
        int O = Order();
        if (O == 0)
            return T();

        // 获得最高次系数，使用累乘方法提高效率
        T Y = m_coefficients.back();
        for (int i = O - 2; i >= 0; i--)
            Y = Y * X + m_coefficients[i];

        return Y;
    }

    // 推入高阶项
    void Add(T k)
    {
        m_coefficients.push_back(k);
    }

    // 多项式降阶，将系数为 0 的高阶项去除
    void Reduce()
    {
        // 只要最后一项非常小，就弹出最后一个元素，除非没有元素
        while (Order() > 0 && fabs(m_coefficients.back()) < 1e-12)
            m_coefficients.pop_back();
    }

    // 多项式积分
    T Intergral(T a, T b) const
    {
        T I = T();
        for (int i = 0; i < Order(); i++)
            I += m_coefficients[i] * (pow(b, i + 1) - pow(a, i + 1)) / (i + 1);
        return I;
    }

    TPolynomial Intergral() const
    {
        TPolynomial poly;
        poly.Add(T());
        for (int i = 0; i < Order(); i++)
            poly.Add(Coeff(i) / (i + 1));
        return poly;
    }

    // 多项式求导
    TPolynomial Derivative() const
    {
        TPolynomial poly;
        for (int i = 1; i < Order(); i++)
            poly.Add(Coeff(i) * i);
        return poly;
    }

    // 获得阶数
    int Order() const
    {
        return m_coefficients.size();
    }

    // 获得系数
    T Coeff(int n) const
    {
        assert(n < m_coefficients.size());
        return m_coefficients[n];
    }
};

//////////////////////////////////////////////////////////////////////////////////////////
///
///

/**
 * 通用多项式类型.
 */
class Polynomial : public TPolynomial<double>
{
  public:
    explicit Polynomial(int o = 0);
    Polynomial(std::initializer_list<double> l);
    Polynomial(const TPolynomial &poly);
    Polynomial(const std::vector<double> &v);
};

//////////////////////////////////////////////////////////////////////////////////////////
///
///

/**
 * 样条类型.
 */
enum class BCType
{
    COMPLETE,   // 完全样条
    NOT_A_KNOT, // 无结样条
    PERIODIC    // 周期样条
};

/**
 * 分段多项式类型.
 */
class PwPolynomial
{
  protected:
    std::deque<double> m_interval; // 分段区间端点
    std::deque<Polynomial> m_poly; // 多项式列表，数量是端点数 +1，包括左右半开区间上的多项式

  public:
    explicit PwPolynomial(Polynomial poly = Polynomial(0)); // 没有区间端点，只有一个多项式
    PwPolynomial(const PwPolynomial &ppoly);                // 拷贝构造
    PwPolynomial(const std::deque<double> &interval);       // 只设置区间，每段区间上是零
    PwPolynomial(const std::deque<double> &interval, const std::deque<Polynomial> &poly); // 设置区间和多项式

    // 多项式赋值
    void operator=(const PwPolynomial &poly);

    // 多项式加减法
    PwPolynomial operator+(const PwPolynomial &poly) const;
    PwPolynomial operator-(const PwPolynomial &poly) const;
    PwPolynomial operator+(const Polynomial &poly) const;
    PwPolynomial operator-(const Polynomial &poly) const;

    // 多项式乘常数
    PwPolynomial operator*(double k) const;
    PwPolynomial operator/(double k) const;

    PwPolynomial operator*(const PwPolynomial &poly) const; // 多项式相乘
    PwPolynomial operator*(const Polynomial &poly) const;   // 多项式相乘
    PwPolynomial operator/(const Polynomial &poly) const;   // 多项式除法

    double Value(double x) const; // 多项式求值

    void AddRight(double t); // 在右端添加端点，增加一个零多项式
    void AddLeft(double t);  // 在左端添加端点，增加一个零多项式

    Polynomial operator[](int i) const; // 获得指定索引的多项式
    Polynomial &operator[](int i);      // 获得指定索引的多项式引用
    Polynomial Poly(double t) const;    // 获得传入点对应的多项式

    int Size() const;       // 获得端点数
    int Order(int i) const; // 获得指定索引的阶数

    PwPolynomial Intergral() const;  // 不定积分
    PwPolynomial Derivative() const; // 导函数

    int Continuity(int i) const; // 检查端点处的连续性
    int Continuity() const;      // 整体连续性

    friend PwPolynomial merge(const PwPolynomial &p1,
                              const PwPolynomial &p2); // 合并两个分段多项式
};

//////////////////////////////////////////////////////////////////////////////////////////
///
/// 辅助函数

/**
 * 多项式整数幂.
 */
Polynomial power(const Polynomial &poly, int n);

/**
 * 产生节点多项式.
 */
Polynomial nodes(const std::vector<double> &c);

/**
 * 已知节点 c0 c1 ... cn ，获得 cj 对应的初等多项式.
 */
Polynomial lagrange(const std::vector<double> &c, int j);

/**
 * (-1,1) 上的 Gauss-Legendre 多项式.
 */
Polynomial legendre(int s);

/**
 * (0,1) 上的首一平移 Gauss-Legendre 多项式.
 */
Polynomial monic_legendre(int s);

/**
 * 获得 Bernstein 基函数.
 */
Polynomial bernstein(int n, int i);

/**
 * 多项式插值.
 */
Polynomial interpolate(const std::vector<double> &x, const std::vector<double> &f);

/**
 * NevilleAitken 算法计算插值点的值.
 */
double neville_aitken(const std::vector<double> &x, const std::vector<double> &f, double t);

/**
 * Hermite 牛顿插值.
 */
Polynomial newton_interp(const std::vector<double> &x, const std::vector<std::vector<double>> &f);

/**
 * 计算差分表.
 */
std::vector<std::vector<double>> calculate_table(const std::vector<double> &x,
                                                 const std::vector<std::vector<double>> &f, int width);

/**
 * 获得 B 样条基函数.
 */
PwPolynomial basic_bspline(const std::deque<double> &knot);

/**
 * 普通 3 次样条插值.
 */
PwPolynomial pp_form_spline(const std::vector<double> &x, const std::vector<std::vector<double>> &f, BCType type);

/**
 * 2 次基数 B 样条在半整数点 1 1.5 2.5 ... 5.5 6 插值，只需要端点值.
 */
PwPolynomial cardinal_bspline_2(const std::vector<double> &x, const std::vector<double> &f);

/**
 * 3 次基数 B 样条在整数点 1 2 3 ... 4 5 6 插值，需要两端一阶导信息.
 */
PwPolynomial cardinal_bspline_3(const std::vector<double> &x, const std::vector<std::vector<double>> &f);

} // namespace CAGD
